Molecular Dynamics Simulations of Temperature Equilibration in Dense Hydrogen 
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The temperature equilibration rate in dense hydrogen (for both Ti > Te and Ti < T^) has 
been calculated with molecular dynamics simulations for temperatures between 10 and 600 eV and 
densities between lO^'^/cc to 10^*/cc. Careful attention has been devoted to convergence of the 
simulations, including the role of semiclassical potentials. We find that for Coulomb logarithms 
£ > 1, a model by Gericke-Murillo-Schlanges (GMS) [Gericke et al., PRE 65, 036418 (2002)] based 
on a T-matrix method and the approach by Brown-Preston- Singleton [Brown et al., Phys. Rep. 410, 
237 (2005)] agrees with the simulation data to within the error bars of the simulation. For smaller 
Coulomb logarithms, the GMS model is consistent with the simulation results. Landau-Spitzer 
models are consistent with the simulation data for £. > 4. 



I. INTRODUCTION 

The strong temperature dependence of thermonuclear 
reaction rates suggests that even small deviations from 
equilibrium can yield differences in burn rates. Thus, the 
pursuit of ignition in the laboratory will benefit from ac- 
curate models of relaxation processes in hot, dense plas- 
mas. One of the greatest uncertainties in the nonequi- 
librium energy balance is the electron-ion temperature 
relaxation rate. Although there have been indirect mea- 
surements for cool dense matter [T], there is no experi- 
mental data in the regime of interest. Even worse, theo- 
retical descriptions of Coulomb collisions suffer from di- 
vergences that make detailed models difhcult to develop. 
Here we take a complementary approach to hot, dense 
plasmas by using molecular dynamics (MD) techniques. 
We use this method to test recent theoretical models and 
compare with standard results. 

The electron-proton coupling rate was first calculated 
by Landau [2] and Spitzer (LS) ^ for classical plasmas 
with weak collisions. They write the electron-proton tem- 
perature exchange rate (l/rpe) in the form, 

Tpe SniempC^ \ vrieC^ mpC^ J Jls ' 

(1) 

where Jls is the LS pre-factor, Ue (up) are the electron 
(ion) number densities, Z = 1 is the proton charge, 
(Tp) are the electron (ion) temperatures, and fcs is the 
Boltzmann constant. C is the so-called Coulomb loga- 
rithm containing details of the collision process. LS used 

Cls = ln(6 max/^min) (2) 

where bmax and bmin are impact parameter cutoffs 
needed to remove divergences that arose from their treat- 
ment, bmin is chosen to be a minimum impact pa- 
rameter consistent with plasma conditions, such as the 
classical distance of closest approach {be = Ze^/ksT). 
At high temperatures, bmin is often modified to include 
quantum diffraction effects by introducing the length 



scale of the electron thermal deBroglie wavelength A = 
y/2TrtP /niekBTe. Typically bmax is chosen to be a screen- 
ing length arising from collective plasma phenomena, 
such as the Debye length Xd = \/ kBTe/AiTe^ne. 

The presence of ad hoc cut-offs and other inconsisten- 
cies led researchers to derive kinetic equations without 
cut-offs [51 El [3 E] • The essence of these theories is the in- 
clusion of strong scattering in the presence of dynamical 
collective (screening) behavior. Two such theories have 
been recently proposed: Gericke, Murillo and Schlanges 
(GMS) ig and Brown, Preston and Singleton (BPS) [9]. 
GMS applied these ideas to dense plasma temperature 
equilibration. They investigated various approximations 
in evaluating of £, including issues with trajectories and 
cutoffs, and provided four different evaluations of the re- 
laxation rate based on quantum kinetic theory. From 
their numerical work, GMS suggest an effective Coulomb 
logarithm [TU] 

Cgmsg = ^ In (1 + [Al, + Ri„] I [AVStt -rbl\), (3) 

1/3 

where Rion = (3/47rnp) ' is the ion sphere radius. This 
expression was described by GMS as the best fit to their 
full T-matrix theory. 

BPS and Brown and Singleton ^\ used dimensional 
continuation to obtain an expression for the electron-ion 
coupling rate accurate to second order in the plasma 
coupling parameter. The method is applicable to both 
degenerate and non-degenerate electrons. For the non- 
degenerate case, they derive 

I^BPS = log(Ai5/A) + (log (167r) - 7 - 1) /2, (4) 

where 7 is the Euler constant. 

The most direct method of studying temperature equi- 
libration in the classical limit is with numerical simula- 
tion; strong, collective scattering at all length scales is 
the forte of MD. Hansen and McDonald (HM) |T2] ex- 
plored temperature equilibration in dense hydrogen us- 
ing MD, comparing their results against a LS model with 
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^■LS = ln(27rAD/A). However, the HM simulations in- 
volved a very small number of particles (TV = 128) with 
presumably large error bars. Here, we expand upon their 
calculations to not only reassess the HM result, but also 
compare with the modern approaches of GMS and BPS. 

II. MOLECULAR DYNAMICS: SIMULATIONS 
AND RESULTS 

MD simulations are applied to two-temperature sys- 
tems of charged particles in a cubic cell with periodic 
boundary conditions. The MD is performed with a fully 
parallel code using a basic leapfrog method [13] with the 
Coulomb interaction evaluated by an Ewald summation 
[131 [Ej. Because the classical Coulomb many-body prob- 
lem is unstable for attractive interactions, we employ 
semi-classical potentials that reduce the Coulomb inter- 
action on short length scales in order to prevent unphys- 
ical, deeply-bound states. We tested several forms of 
the diffractive [T71 US] and Pauli [TO] HU] terms for these 
potentials. The resulting equilibration times typically 
vary by less than 15%, which is within the statistical er- 
ror of the MD data. The similarity is not unexpected, 
since most semi-classical potentials resemble one another 
above 10 eV f2?. We report results using the semi- 
classical potential in HM [T^ . 

Vab{r) = ^^(l-exp(-27rr/A,b)) 

+ fcBTln2cxp(-47rrVA^b/ln2)^,e<56e (5) 

where Kab — \/'^t^K^ / HabknT , fiab is the reduced mass, 
and T — T(. except when a and b are both protons when 
T = Tp. The potentials are temperature-dependent, but 
were held constant in most of our short simulations. For 
long simulations to equilibration, wc allow the tempera- 
ture parameters to evolve with time, using a smoothed 
exponential average of the instantaneous MD value. 

Simulations were run long enough to extract a re- 
laxation time (typically 10% of r*), with some strong- 
coupling cases continued to complete equilibration. We 
obtain Tpe by fitting the temperature over a brief interval, 

Tp X'g ^ dTp T^g Tp 

dt ^ep dt ^pe 

We choose the timestep to conserve total energy over 
the entire simulation (^AE/E < 10^'*) when using fixed 
potentials. Typically, At ranges from 5 x 10~^ to 10"'^ fs. 
Any drift in the energy is tightly controlled, as artificial 
heating can distort the true relaxation rate. In long runs, 
the potentials change slowly as the temperature relaxes. 
Although energy is not conserved in these cases, the total 
energy change remains less than 3%. In practice, Tpe 
calculated from the time-dependent potential is within 
10-15% of the result for the constant potential. 

Convergence with respect to system size is tested by 
employing various particles numbers N ranging from 



Case 


ni(l/cc) 


T^eV) Ti{eV) 






A 




10.0 


20.0 


2.04 X 10* 


4.9 X 10^ 


B 


10^° 


30.0 


60.0 


7.89 X 10* 


4.3 X 10* 


C 




10.0 


20.0 


5.23 X 10^ 


1.7 X 10^ 


D 


^q22 


30.0 


60.0 


1.73 X 10^ 


6.6 X 10^ 


E 


^q22 


100.0 


200.0 


6.45 X 10^ 


2.2 X 10^ 


F 




10.0 


20.0 


8.87 X 10^ 


3.5 X 10^ 


G 




30.0 


60.0 


8.27 X 10^ 


3.3 X 10^ 


H 




100.0 


200.0 


1.72 X 10^ 


6.2 X 10^ 


I 




300.0 


600.0 


4.17 X 10^ 


8.0 X 10^ 


J 


L61 X 10^" 


29.9 


80.1 


20.2 X 10^ 


5.3 


K 


L61 X 10^'' 


91.47 


12.1 


1.20 X 10^ 


1.7 X 10^ 


L 




100.0 


200.0 


3.65 X 10^ 


3.2 X lO'^ 


Ml 


^q20 


10.0 


40.0 


2.05 X 10* 


3.0 X 10^ 


Ma 


10^" 


10.0 


40.0 


2.18 X 10* 


4.5 X 10^ 


Mg 




10.0 


40.0 


2.28 X 10* 


9.6 X 10^ 



TABLE I: Density, initial electron, and ion temperature, re- 
laxation time and standard deviation of the MD simulations. 



N = 128 (the number that HM employed), to as many 
as = 64, 000. The results reported here use N = 1024. 
Statistical uncertainty for each case is estimated by com- 
puting the relaxation rate from equivalent samples (from 
8 to 64) of a microcanonical ensemble and then taking 
the average and standard deviation. Sensitivity to initial 
conditions is studied using the ensemble of simulations 
and/or by discarding a portion of the initial temperature 
evolution. 

The nonequilibrium system is prepared using two sep- 
arate Langevin thermostats for protons and electrons. 
Initial configurations are sampled from a stationary dis- 
tribution obtained after 10^ — 10^ timesteps. The ther- 
mostats are then removed, and the species allowed to un- 
dergo (microcanonical) coUisional relaxation for approx- 
imately 10^ timesteps. 

Equations [6] are valid for an ideal gas equation of state 
for the plasma. For strongly coupled plasmas there is a 
significant potential energy contribution that would in- 
validate this assumption. However, the error associated 
with using the temperature evolution equations is small 
in the temperature-density regimes of interest here[T^. 
Although the MD temperature relaxation is asymmet- 
ric in the strong-coupling cases, we find \dT^/dt\ and 
\dTp/dt\ differ by only about 10%. Thus, we only report 

1/r* = 1/Tpe + l/Tep « 2/Tpe. 

Table 1 lists the set of initial conditions for 15 dif- 
ferent systems. The ensemble average temperature re- 
laxation, T* (calculated from dAT/dt = AT/r*), and 
the standard deviation, tr, are in femtoseconds. A range 
of initial conditions were chosen to span the weakly- to 
strongly-coupled and the degenerate to non-degenerate 
regimes. We include two sets of initial conditions con- 
sidered by HM (Cases J and K). In most cases, hydro- 
gen plasma is simulated using the true electron-proton 
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mass ratio of 1:1836. In Case L, the cold electrons were 
replaced with cold protons in order to shorten the re- 
quired simulation time. Cases Mi_3 involve a compari- 
son of electron-proton and positron-proton systems and 
will be discussed below. Cases F and G have degenerate 
electrons. Degeneracy effects are treated in neither the 
classical MD simulations nor in the LS, GMS6, or BPS 
models. Hence, the models can be directly compared to 
the simulations even for those cases when comparisons 
with experiment would be questionable. 



III. COMPARISON WITH THEORY 

Figure [T] shows MD results for Case K run to near- 
full relaxation using potentials that are implicitly time- 
dependent (temperature-dependent). We also display 
predictions for LS, GMS6, and BPS. The MD data in 
Figure [l] is most closely matched by BPS (although this 
is partly fortuitous, as will been seen below) followed by 
GMS6. The LS model predicts the fastest relaxation, ex- 
ceeding MD by about a factor of two. This disagreement 
contradicts the conclusion reached by HM. At the same 
time, our r* for Cases J and K agree with those reported 
by HM. We attribute the discrepancy to inconsistent def- 
initions of Tpe, and T*: ti^s is properly equal to Tpe, 
which is 2r* (not r*) if Tpe = Tep. As previously noted 
by HM, however, ambiguities in the 6mm and bmax may 
be sufficient to accommodate this difference. 

To make comparisons of our MD results with theoret- 
ical predictions more transparent, we define an effective 
Coulomb logarithm as Cmd = "^JlsIt*- This result is 
then compared with the theoretical prediction for C com- 
ing from LS, GMS6, and BPS. Fig. Q shows simulation 
results for Cmd with error bars along with theoretical 
predictions for Cgmsq (solid) and Cbps (dashed) as a 
function of initial electron temperature. Numerical re- 
sults and analytic expressions for C are arranged accord- 
ing to density; n = 10^°, 10^^, and 10^^ (blue, red and 
black respectively). 

In regions where it is expected to be applicable, we 
find that LS systematically overestimates the effective 
Coulomb logarithm and thus predicts a relaxation rate 
that is too fast relative to the MD results. For plasmas 
with £ > 1, the MD results are consistent with both 
the GMS6 and BPS, suggesting that approaches beyond 
LS are indeed more predictive. As expected, BPS in- 
creasingly underestimates the relaxation rate for £ < 1; 
BPS is not intended for use in this regime. For the case 
shown in Fig. [T] the underestimation at lower tempera- 
tures compensates for an overestimation at early times, 
making agreement with this simulation fortuitously good. 
As is evident from Figure [2] this would not be the case in 
general [23] . We find that GMS6 captures the qualitative 
variation of £ over a surprisingly broad range of density 
and temperature. Further discrimination between these 
theories in the region where where they are expected to 
be most accurate (low density and high temperature) is 
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FIG. 1: Electron (top curves) and proton (bottom curves) 
temperature relaxation is shown based on MD, GMS, LS, and 
BPS for Case L. The MD results are shown by points from 
several simulations, with a line through the average. Note 
that all approaches relax slower than LS. 




Temperature (eV) 



FIG. 2: Theoretical (GMS6 [solid], BPS [dashed], and LS 
[dotted]) and MD calculations of £ as a function of initial 
Te for densities 10^"/cc, lO^Vcc, and 10^''/cc (blue, red and 
black respectively). Additional detail is in the text. 
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not possible given the large uncertainties present in our 
current MD simulations. However, our results suggest 
that validation of these theories could be accomplished 
with carefully controlled experiments [24^ and larger (and 
longer) simulations that further reduce statistical error. 

Finally, LS predicts identical equilibration rates for 
like-charge and opposite charge systems. We tested this 
by performing three sets of simulations at the same den- 
sity and temperatures (Case Mi_3 in Table I.) We simu- 
lated electrons-protons (Mi) and positrons-protons (M2) 
using Equation 6, and positrons-protons using a pure 
1/r Coulomb potential (M3). The relaxation rates for 
all three cases agree to within our error bars, suggest- 
ing that energy transfer in these systems is occurring 
predominately on length scales longer than the thermal 
deBroglie wavelength. 

IV. CONCLUSIONS 

We have performed MD simulations of the tempera- 
ture relaxation process in hot, dense hydrogen. We in- 
vestigated systems containing as large as 64,000 parti- 
cles, finding that N ~1000 particles is sufficient for most 
cases we considered. Our simulations span a large range 
of temperature and density parameter space, including 
the first simulations in the low-density, high-temperature 
limit. 



For the weakly coupled plasmas where C>1, the sim- 
ulations are consistent with both GMS6 and BPS. In con- 
trast, the LS approach systematically overestimates the 
relaxation rate. In the limit of high temperature and low 
density, all models arc in agreement, however. Our MD 
results suggest that LS is accurate for £ > 4, rather than 
the usual restriction oi C > 10, in agreement with pre- 
vious work [8, 21 . More modern approaches exemplified 
here by GMS6 and BPS clearly extend the accessible pa- 
rameter space closer to C ~ 1, with GMS6 providing a 
reasonable description of the MD data even for for £ < 1. 

We have employed two forms of the semiclassical po- 
tentials needed for stability in an MD simulation with 
attractive potentials, and have found a very slight effect 
from the form of the potential; as such, we believe that 
our results are not sensitive to the choice of the semiclas- 
sical portion of the potentials. 
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